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We calculate ground-state properties of a many-quark system in the string-flip 
model using variational Monte Carlo methods. The many-body potential energy of 
the system is determined by finding the optimal grouping of quarks into hadrons. 
This (optimal) assignment problem is solved by using the stochastic optimization 
technique of simulated annealing. Results are presented for the energy and length- 
scale for confinement as a function of density. These results show how quarks cluster- 
ing decreases with density and characterize the nuclear- to quark-matter transition. 
We compare our results to a previously published work with a similar model which 
uses, instead, a pairing approach to the optimization problem. 
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I. INTRODUCTION 



Hadronic models have been very successful in describing the interactions which give rise 
to the properties of nuclei and nuclear matter. However, at a more fundamental level, 
our current understanding of the strong force gives us a picture of fermion fields (quarks) 
interacting via gauge fields (gluons) generated by an underlying local SU(3) C symmetry. 
Perhaps the greatest challenge in the field of nuclear structure is to understand the emergence 
of nuclear properties from the underlying QCD dynamics in which quark and gluon fields 
are the fundamental degrees of freedom. The success of experimental efforts at CEBAF and 
RHIC aimed at identifying signatures of the quark substructure in nuclei will rely on the 
ability of theoretical models to identify the appropriate observables. One might be interested, 
for example, in understanding how are hadronic properties such as the nucleon form factor 
modified by its interactions with the medium? Are there new modes of excitation in many- 
quark systems (e.g., coherent oscillations in density, spin, flavor, or color) not present in 
purely hadronic models? Also, what is the nature of the nuclear-matter to quark-matter 
transition? The hope is that simple models such as the one to be discussed here will provide 
a first step towards developing "realistic" many-quark models which identify the relevant 
nuclear observables. 

In this paper we will be working with a string-flip model which has been described 
in detail elsewhere ■ For the many-quark system, the model attempts to reproduce 
nuclear /quark- matter properties using exclusively quark degrees of freedom. Confining 
strings which adjust instantaneously to the changing constituent quark positions (adiabatic 
approximation) are used as a minimum representation of the extremely complex quark-gluon 
dynamics. In this model, there are no strings connecting quarks "belonging" to different 
hadrons. Hence, the force saturates, thus allowing clusters to separate, and avoids the emer- 
gence of long-range van der Waals forces. Among the desirable properties of this potential 
are explicit quark-exchange symmetry, asymptotic freedom and confinement. However, be- 
cause the model is a many-body generalization of the nonrelativistic constituent quark model 
it is not chirally symmetric, explicitly gauge invariant, nor is it Lorentz invariant. 

A critical feature shared by all models of this type is the need to determine an optimal 
grouping of the quarks which minimizes the total potential energy stored in the strings. 
This assignment problem is similar in nature to the well known traveling salesman problem 
in the field of combinatorial optimization. Hence, except for simulations with a very small 
number of quarks, the computational cost of an exhaustive search over the possible string 
configurations is beyond the capabilities of even the most powerful computers available today. 
As a result, we present a model which uses the technique of simulated annealing to solve the 
optimal grouping problem. For the purposes of comparison, we will also describe a similar 
model which uses an efficient pairing algorithm to solve a simpler assignment problem ||. 

The paper is organized as follows. Section || gives a general description of the string-flip 
model. In this section the many-body potential is introduced as well as the techniques used 
in solving the optimization problem. In Sec. [Hj we define the one-parameter variational 
wavefunction that is employed to generate, via Metropolis Monte Carlo, the quark config- 
urations. Results are presented in Sec. [IV| while Sec. [V] summarizes the main results and 
points toward possible future work in this 
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II. STRING-FLIP MODEL 



A. The Many-Body Potential 

In the present version of the string-flip model the quarks are grouped into colorless 
triplets (red+green+blue=white) by harmonic confining strings [|J . Governing the dynamics 
of N = 3A quarks grouped into A colorless hadrons is the following Hamiltonian, 

N V 2 

H = V--^+V, (1) 

1=1 



where 



V = Min£ |[(xf - xf) 2 + (xf - xf) 2 + (xf - xf f\ . (2) 

•t=i 

Here xf , xf , and xf are the positions of the red, green, and blue quarks belonging to hadron 
i and the minimization is taken over a (potentially enormous) set of {A\) 2 possible string 
configurations. Furthermore, since the quark mass (M) and the spring constant (k) are the 
only dimensionfull quantities in the problem one is free to set 

k = M = 1 , (3) 



and, thus, measure all energies in units of the oscillator frequency Jk/M, and all lengths 

according to the oscillator length (Mk)~^. It is then easy to rescale any results according 
to a particular choice of k and M. 

The potential energy for an isolated 3-quark cluster (nucleon) is, in principle, an arbitrary 
confining function of the three quark coordinates. In the simplest version of the model used 
here, however, quarks are assumed to be confined by harmonic strings. For this case a 
triangular (or A) configuration yields identical results to the Y configuration which has 
quarks connected to a three-quark junction as seems to be preferred by the underlying 
gauge symmetry of QCD (see Fig. |I|). We have chosen the A configuration since it simplifies 
the evaluation of the potential energy by eliminating the calculation of a center-of-mass 
coordinate for each hadron. 

In this model there are no long range van der Waals forces between hadrons. The only 
(residual) interactions among hadrons are generated by the possibility of quark exchange 
(change in grouping) and the Pauli exclusion principle between quarks. Thus, the residual 
interaction between hadrons is short range (of the order of the root-mean-square radius of 
the hadron). Furthermore, the potential is a truly many-body operator (i.e., it cannot be 
reduced to a sum over two-body operators); moving a single quark may cause all the quark 
groupings to change. 



B. The Three-quark Assignment Problem 

For the many-quark system, local SU(3) C gauge symmetry demands that a color flux 
tube (string) leaving one quark must connect to an anti-quark or to a three-quark junction. 
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For harmonic strings one can avoid the use of the Y-coupling scheme and can connect the 
quarks directly to one another. So, in our model, a string leaving one quark will connect to 
another. The fundamental question is which one? Lattice QCD solves this assignment prob- 
lem automatically by evolving gluonic degrees of freedom according to extremely complex 
dynamics resulting in various flux tube arrangements. Unfortunately, the computational 
cost of solving full QCD on a lattice is enormous. Perhaps for the forseeable future approx- 
imate models of the complex QCD dynamics will be needed to give insight into the role of 
hadronic substructure in nuclear physics. In the string-flip model used here the requirement 
is to determine the string topology which minimizes the potential energy. The flux tubes 
are assumed to adjust instantaneously to the changing constituent quark positions (adia- 
batic approximation) and are used as a minimum representation of the extremely complex 
quark-gluon dynamics. 

For the problem of grouping N = 3A quarks into A colorless hadrons, an exhaustive 
search for the minimal configuration has a computational cost proportional to A\ 2 . This 
becomes prohibitive for more than about 5 hadrons (note, for ten hadrons there are 10 13 
possible groupings!). Another approach is obviously needed. 

1. Simulated Annealing 

To date, no efficient (power-law) algorithm has been developed to solve the three-quark 
assignment problem. Thus, in order to evaluate the many-quark potential [Eq. (0)] we resort 
to simulated annealing. 

Simulated annealing is a stochastic optimization technique which has been successfully 
applied to problems of large dimension. For background on the annealing approach to 
problems of combinatorial optimization, we refer the reader to the excellent book by Otten 
and van Ginneken ||. In our problem the quantity to be minimized is the potential energy 
and we choose the search space to be all colorless triplets of quarks. Here, iV = 3A, and 
there are A\ 2 possible string configurations to search over. 

The annealing simulation proceeds much like the familiar Metropolis random walk used in 
numerical integration except for the addition of a control parameter, or temperature, which 
is initially set very high. A high temperature is one which is much greater than typical 
energy differences encountered by the random walker when taking steps. This means that 
the probability for accepting the next step is nearly unity and the walker wanders randomly 
throughout the search space, uninhibited by the structure of the energy landscape. The 
control parameter is then slowly (adiabatically) lowered so that the walker progressively 
visits lower energy regions with increasing probability. By 'slowly' we mean that the system 
must always remain near equilibrium (quasi-equilibrium). Eventually, as the temperature is 
taken lower and lower, the walker 'freezes' into what is hopefully the lowest possible energy 
configuration. The probability that the random walker 'freezes' into this global minimum 
instead of one of the innumerable local minima is directly related to how well the criterion of 
quasi-equilibrium was maintained throughout the cooling process and this in turn is related 
to the amount of computation time invested. So, in practice, one always solves the problem 
with some confidence level determined by computational constraints. For example, in the 
string-flip model, the estimate of ground-state observables will be done via Metropolis Monte 
Carlo. Thus, the annealing must be performed at each step of the Monte Carlo simulation 
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of the quark coordinates. This means an enormous number of annealing runs must be done 
during a typical simulation. Obviously there will be limits on how long one can spend on 
each annealing process. 

The simple annealing algorithm used in this work is described in Fig. |2|. By far the most 
important issue in developing this kind of simulation is the selection of a cooling schedule 
and choosing the number of thermalization steps. In general, the number of thermalization 
steps could be a function of temperature. In this particular simulation, we have chosen a 
constant number of thermalization steps. The initial temperature should be chosen high 
enough so that the entire configuration space is sampled but not so high that computation 
time is wasted. Also, the cooling must be slow enough that the random walker doesn't 
prematurely freeze into a local energy minimum but not so slow that computation time is 
unnecessarily long. 

Implementation of the annealing algorithm within the string-flip model is complicated 
by the fact that it is imbedded in a larger Monte Carlo simulation (i.e., calculation of the 
expectation values (V)) of the quark coordinates. Each time the configuration of quarks is 
modified, the energy landscape for the subsequent annealing process changes slightly and it 
is the structure of the energy landscape which drives the choice of cooling parameters. The 
energy landscape will also vary significantly as a function of density. For example, in the low 
density hadronic phase where the quarks are tightly clustered, there is a large energy contrast 
between states in the vicinity of the global minimum and states which lie farther away in 
the configuration space. Now, as the density increases, the system experiences a transition 
to a phase characterized by somewhat larger quark clusters. At these higher densities, the 
energy difference between random groupings and those near the global minimum is less than 
in the low density phase. We can see then that a cooling schedule should be robust enough 
to perform adequately at all densities of interest. We chose our cooling parameters based on 
simulations in the high density region where the annealing is somewhat more demanding. 

Inevitably, a lot of numerical experimentation is required to determine a good set of 
cooling parameters for a new problem. In the work presented here, we have chosen an 
exponential cooling schedule and a number of thermalization steps that yields a confidence 
level of 95 per-cent for finding the global minimum. This confidence level was determined 
by compiling statistics from many simulations involving 18 randomly distributed quarks. In 
each case, we compared configurations arrived at by an exhaustive search with those given 
by the annealing algorithm to see if the global solution was found. 

Figure gives an example of how the simulated annealing simulation progresses on 
average for N = 24 and a low density p = 0.05. This is in the hadronic phase of the system. 
The average is over many Monte Carlo generated instances of the quark coordinates (quark 
coordinates are generated according to a variational wavefunction that will be described in 
Sec. HI). The temperature, acceptance ratio, average potential energy, and best potential 
thus far are plotted as a function of the cooling step. A typical plot of an individual annealing 
process would show an abrupt drop in the energy (freezing) somewhere between steps 40 
and 55, a common feature of simulated annealing runs. 
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2. Pairing Algorithms 



One way to get around the grouping problem described above is to consider an alter- 
native, albeit simpler, grouping scheme. In this scheme one considers, instead, pairing the 
quarks two colors at a time. Economists have long been interested in the problem of pairing 
N factories with N retail stores in order to minimize the overall cost of exchanging goods. 
Efficient pairing algorithms with a computational cost proportional to N 3 have already 
been developed by mathematicians 0. Adapting these efficient algorithms to the problem 
of interest implies calling the pairing algorithm three times in succession in order to indepen- 
dently pair red and green quarks, green and blue quarks, and blue and red quarks. Notice, 
however, that the many-quark dynamics must be modified in order to take advantage of the 
pairing algorithm. The potential energy for the system is now given by, 

V = Min^ \{xf - x?Y + Mm^-ixf - xff + Min^ \{xf - xff , (4) 
i=i A i=i A i=i A 

and should be compared to Eq. (^). This version of the string-flip model has already been 
used in Ref. @] and we include some results here for comparison. 

In spite of this change both models share many common features. Clearly, they are 
identical in the case of an isolated cluster. Moreover, both guarantee cluster separability 
thus avoiding the emergence of long-range van der Waals forces. There are, however, some 
differences. Most notoriously, by pairing quarks independently there is no guarantee that 
the quarks will be grouped into three-quark clusters. Color-neutral hadrons in this model 
may contain any multiple of three quarks. Thus, the space of string configurations over 
which the potential is minimized consists of all colorless clusters containing multiples of 
three quarks. We refer to this model as the multi-quark-cluster (MQC) model. In contrast, 
the model based on the potential given by Eq. (Q), which allows for only three-quark clusters, 
will be referred as the three-quark-cluster (TQC) model. Since the minimization space is 
larger in the MQC model than if only 3-quark clusters were allowed, one can generally 
expect solutions which are lower in energy. However, this will not be the case in the low- 
density regime, where hadrons are well separated and pairing quarks independently should, 
nonetheless, generate only three-quark hadrons. 



III. THE VARIATIONAL WAVEFUNCTION 

In this work we use a simple single-parameter variational wavefunction to represent the 
ground state of the many-quark system, 

^ = e~ AV * FG . (5) 

The Fermi-gas wavefunction, ^fg^ is exact for a system of free fermions with no correlations 
other than those generated by the Pauli exclusion principle. It is a product of red, green, and 
blue Slater determinants. The exponential factor expresses the degree of clustering through 
the variational parameter A with A" 1 / 2 being the length-scale for quark confinement. In 
spite of the simplicity of the wavefunction, it is exact, both, in the low-density nuclear 
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matter phase where the isolated-hadron limit is recovered (A = 1/V3) and in the high- 
density Fermi-gas phase where A = 0. The wavefunction regards red, green, and blue quarks 
as distinguishable particles (in the same sense as spin-up and spin-down electrons can be 
regarded as distinguishable in the presence of a spin-independent Coulomb interaction). 
Thus the Pauli principle is enforced only among quarks of the same color. 



A. Determination of A 

The ground state expectation values for the kinetic and potential terms of our Hamilto- 
nian can be related using integration by parts. The result is, 

(tf|T|tf) = T FG + 3A 2 (^|V|^) . (6) 

Tfg is the kinetic energy of a free Fermi gas. The term 3A 2 (V) is the increase in kinetic 
energy above the Fermi-gas value generated by clustering correlations. Thus, to calculate 
the total energy we only need the expectation value of the potential, 

(#|H|tf) = E(X) = T FG + (3A 2 + |V|tf) . (7) 

As in any variational approach, the variational parameter is determined by minimizing 
the energy. A 'brute force' implementation is to simply calculate E(X) for several values of 
A and use this information to estimate A m i n . This is very expensive computationally since 
several Monte Carlo estimates of E(X) are required for each estimate of A min . 

Fortunately, there is a more elegant and efficient method for extracting A m i n which takes 
maximum advantage of every Monte Carlo run. It makes use of scaling relations for the 
potential when it is evaluated at two different values of A and the quark density p. The 
condition for finding A min is the usual vanishing of the derivative of the energy with respect 
to the variational parameter 

dE 

— = 6A(V)-[6A 2 + 2][(V 2 )-(V) 2 ]=0, (8) 

which requires, in addition to the potential, an estimate of the expectation value of the 
square of the potential. 

At first glance, (V) appears to be a function of both A and p. However, on the grounds 
of simple dimensional analysis one can show that the potential energy can be written (for a 
fixed number of quarks) as 

(V) Ap /L 2 = v(XL') , (9) 

where p = N/L 3 and v is a dimensionless function of the scaling variable AL 2 . Hence, we see 
that (y)\ p /L 2 is in reality a function of only one scaling variable leading to the following 
scaling relation for the potential: 

(V)a p =(^(V) av , (10) 

with 
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A' = (-) f A. (11) 
P 

We can use this freedom as follows. After performing a calculation using p, A we can 
always rescale our results to a new pair of variables, p', A', such that dE/dX = is satisfied. 
After some algebra, we find the new density which accomplishes this is given by, 



^ i3 v>r-w Ari - <i2) 



This is indeed a useful result since it now takes only one Monte Carlo estimate of (V) and 
(V 2 ) to determine the variational parameter. The only potential pitfall of this method for 
acquiring A m ; n is that E(X) should be a smooth function without multiple minima. One 
should therefore be careful in using this technique, particularly near the transition between 
the confined and unconfined phases of the system, where double minima or degeneracies 
might occur. 



IV. SIMULATION RESULTS 

Figure |] shows the difference in the potential energies in the MQC model (obtained by 
using the pairing algorithm) and the TQC model (with simulated annealing) as a function of 
the density of the system. As expected, the two agree at very low density where the hadrons 
are well separated and multi-quark configurations are suppressed. At a higher density, 
however, the MQC model takes advantage of the possibility of forming large multi-quark 
rings and yields a lower energy. 

Figure shows the energy per quark as a function of density for the MQC model (N = 
24, 96) and for the TQC model (N = 24) along with the Hartree-Fock (A = 0) result. 
Deviations from the Hartree-Fock result in the high-density phase are indicative of finite- 
size effects present in the simulation and get smaller as N is increased in the MQC model 
as one might expect. While this finite-size effect lowers the energy per quark in the quark- 
matter phase, it has little effect in the nuclear-matter phase. The transition between the 
two phases occurs abruptly near p = 0.12. 

Figure [| also reveals that the TQC model with N = 24 shows roughly the same finite-size 
effect as the N = 96 MQC result. One might expect larger finite-size effects for the latter 
since the MQC model allows the formation of large multi-quark clusters which can have 
dimensions comparable to the size of the box. In contrast, we have found that three-quark 
clusters are always considerable smaller than the size of the box. Based on this evidence, 
we will compare N = 24 annealing results with iV = 96 pairing results since they seem to 
have comparable finite size effects. 

The variational parameter A as a function of density is given in Figs. |6], [7], |8|. These 
figures show the transition between the hadronic phase where A is near the isolated hadron 
value, and the quark-gas phase characterized by A tending to zero. The MQC model shows a 
step-like behavior in X(p) at the transition density, whereas the TQC model gives a gradual, 
continuous transition. This can be understood in terms of the structure of E(X) near the 
phase transition. E(X) generated by the MQC model is characterized by the formation 
of double minima near the transition giving rise to an abrupt shift to a smaller value of 
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A . In contrast, we have found that in the annealing model, E(X) becomes very flat 
(constant) near the phase transition, resulting in large fluctuations in the string lengths and 
a corresponding broad and smooth transition between the two phases (the nature of the 
transition for a larger number of quarks in the TQC model remains to be investigated). 

The phase transition is demonstrated most dramatically in Fig. ^| where the abrupt 
increase in the potential energy (increasing string length) and a rapid decrease in kinetic 
energy (lower Fermi energy, deconfinement) can be seen at the crossover point. At very low 
density the equipartition of energy between kinetic and potential forms is satisfied, and at 
high density, the system is dominated by the Fermi-driven kinetic component as expected 
from the behavior of a free quark gas. In Figure |1(] the average string length was directly 
measured as a function of density. The length is presented in units of the isolated cluster 
value which can be calculated exactly. Again, we see the abrupt swelling of the clusters near 
the transition density. As expected, the behavior of the average string length is identical to 
that of the potential energy in Fig. [| 

V. CONCLUSION 

In this paper, we have presented the main results of a three-quark string-flip model using 
simulated annealing to determine optimal groupings. Previous work with a similar model 
used pairing to achieve the grouping. A major difference between the models is the presence 
of multi-quark clusters in the MQC model, while the TQC model allows for only 3-quark 
hadrons. 

To summarize, both models have the expected low density (isolated hadrons) and high 
density (quark-gas) behavior. Also, the transition density is roughly the same in both 
models. The main differences are that the TQC model appears to have smaller finite size 
effects due to the 3-quark constraint, and the variational parameter, A(p), is smooth and 
continuous across the phase transition whereas, in the MQC model, A(p) is step-like and 
discontinuous. 

The immediate followup to this work will be to use the values determined here for the vari- 
ational parameter, A, to calculate additional ground-state observables like the quark-quark 
correlation function. This might further characterize the hadron to quark-gas transition 
and could be important in calculating 'quark giant resonances'. Longer range projects may 
include a spin-dependent interaction (which is responsible for the N — A mass splitting) and 
light-quark flavor degeneracy in the hope of understanding the saturation of nuclear matter 
using only quark degrees of freedom. In addition, one might include additional (heavy) 
flavors to study and characterize the transition to strange matter. 
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FIGURES 



FIG. 1. A and Y coupling schemes for quarks. For harmonic potentials, they are equivalent. 

FIG. 2. Annealing algorithm used for solving the string assignment problem. 

FIG. 3. Various quantities involved in the annealing process as a function of the number of 
cooling steps. This figure is an average result over many Metropolis steps for 24 quarks. 

FIG. 4. Percent difference between optimal potential energies as returned by the TQC model 
(annealing algorithm) and the MQC model (pairing algorithm). The larger configuration space of 
the pairing process yields lower energies except at very low density where the difference tends to 
zero. 

FIG. 5. Comparison of the total energy per quark in the TQC model vs. the MQC model for 
24 and 96 quarks. Also shown is the Hartree-Fock (A = 0) result. 

FIG. 6. Variational parameter vs. density from the MQC model with N = 24 quarks. 

FIG. 7. Variational parameter vs. density from the MQC model with N = 96 quarks. 

FIG. 8. Variational parameter vs. density from the TQC model with N = 24 quarks. 

FIG. 9. Results of the TQC model using simulated annealing for 24 quarks. Total, kinetic and 
potential energy per quark are shown together with the Hartree-Fock (A = 0) result. The vertical 
dashed line indicates the transition region between the hadronic and quark phases. 

FIG. 10. Average string length (units of isolated hadron value) as a function of density in the 
TQC model. 
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